%pylab inline
import matplotlib.pyplot as plt

from scipy.integrate import * # az integráló rutinok betöltése

import warnings
warnings.filterwarnings('ignore')

Populating the interactive namespace from numpy and matplotlib

Debye model for the specific heat capacity:

$C_V= 9 p N k_{\mathrm{B}} {\left(\frac{T}{\Theta_\mathrm{D}}\right)}^3 \int_0^{\Theta_\mathrm{D}/T}\, \frac{x^4 e^x}{{\left(e^x-1\right)}^2}\, dx$,

where $p$ is the number of atoms in the basis and $N$ is the number of unit cell in the sample.

For low temperature $C_V \approx p N k_{\mathrm{B}}\, \frac{12 \pi^4}{5} {\left(\frac{T}{\Theta_{\mathrm{D}}}\right)}^3$,

when $T\ll \Theta_\mathrm{D}$.

See Eq. (12.3.35) in Jenő Sólyom: Fundamentals of the Physics of Solids, Volume 1: Structure and Dynamics (A modern szilárdtestfizika alpajai I. A szilárd testek szerkezete és dinamikája) see here on page 416.

# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 9   #    12 nagy abrahoz
yfig_meret= 6    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 20  #  30 nagy abrahoz
legend_meret= 20   #  30 nagy abrahoz

def ct(t):
    expt= exp(t)
    res = t**4*expt/(expt-1)**2
    return res

Np = 100

Tmax=1.7
pici = 0.001

T = linspace(pici,Tmax,Np)
Tm = linspace(pici,0.17,Np/2)
T3 = 12*pi**4/5*Tm**3

cV=[]
for tt in T:
    cV.append(9*tt**3*quad(ct,pici,1/tt)[0])


figsize(xfig_meret,yfig_meret)
plot (T, cV,color='r')
plot (Tm,T3,
      label=r'$\frac{C_V}{p N k_{\mathrm{B}}} \approx \frac{12 \pi^4}{5} {\left(\frac{T}{\Theta_{\mathrm{D}}}\right)}^3$')


xlabel(r'$T/\Theta_{\mathrm{D}}$',fontsize=xylabel_meret)
ylabel(r'$C_V/(p N k_\mathrm{B})$',fontsize=xylabel_meret)

legend(loc='center right',fontsize=legend_meret)

title(r'$\mathrm{Debye \,\, fajhő}$',fontsize=20)


xlim(0,Tmax)

grid();

#savefig('Fig_Debye_fajho.pdf');  # Abra kimentese